#######		REPLICATION FILES
#######		Climate Variability and Irregular Migration to the European Union
#######		Global Environmental Change
#######		Fabien Cottier and Idean Salehyan
#######		Replication: Figure 2: map displaying number of irregular migrants 2010-2015
#######		This version: April 2021




# Set working directory
setwd("path/to/replication/directory")


# Load packages
library(sf)
library(dplyr)
library(tmap)
library(raster)
library(tidyverse)
library(sp)
library(rgdal)
library(WDI)


# specify destination projection: Robinson projection
proj="+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
proj_name="_rob"


# load cshapes country borders dataset and reproject
cshapes_proj <- sf::read_sf("shp/cshapes_mollenweide.shp")
cshapes_proj <- cshapes_proj[cshapes_proj$COWEYEAR>=2015 & cshapes_proj$COWSYEAR<=2015,]
cshapes_proj <- sf::as_Spatial(cshapes_proj)  # cumbersome, but only way to correctly reproject maps without errors...
cshapes_proj <- sp::spTransform(cshapes_proj, proj)
cshapes_proj <- st_as_sf(cshapes_proj)


### add frontex data and create indicators for UK, Ireland + Balkans countries
frontex <- readr::read_csv("FrontexIBCP_2009_2017m_cntr_df.csv")
frontex <- frontex %>%
  dplyr::filter(year %in% 2010:2015) %>%
  dplyr::group_by(cowc) %>%
  dplyr::summarise(
    nationality = first(nationality),
    nmigr = sum(nmigr),
    nmigr_exbalk = sum(nmigr_exbalk),
    frontex_data = first(frontex_data)
  )
cshapes_proj <- dplyr::left_join(cshapes_proj, frontex,by=c("COWCODE"="cowc"))
cshapes_proj <- cshapes_proj %>% 
  dplyr::mutate(
    non_eu=ifelse((!is.na(nmigr) & is.na(nmigr_exbalk)) | (COWCODE %in% c(200, 205)), 1, 0)
  )
cshapes_proj_migr <-  cshapes_proj[!is.na(cshapes_proj$nmigr_exbalk),]
cshapes_proj_non_eu <- cshapes_proj[cshapes_proj$non_eu==1,]
cshapes_proj_euefta <- cshapes_proj[cshapes_proj$non_eu==0 & is.na(cshapes_proj$nmigr_exbalk),]


### create graticule for lon / lat
grid <- sp::gridlines(
    Spatial(
        matrix(c(-180, 180, -90, 90), 2, 2,byrow = T, dimnames = list(NULL, (c("min", "max")))),
        CRS("+init=epsg:4326")
    ),
    easts =  seq(-180, 180, 45),
    norths = seq(-90, 90, 30),
    ndiscr = 1000)
grid_proj = sp::spTransform(grid, proj)


## create map
legend_N<-"N irregular migrants"
map_N <- tm_shape(grid_proj) +
    tm_lines(col = "grey", alpha = .3) +
    tm_shape(cshapes_proj_migr) +
    tm_polygons(
        col = "nmigr_exbalk",
        breaks=c(0,10^c(2:6)),
        lwd=.35,
        title=legend_N,
        legend.is.portrait=F,
        interval.closure = "left",
        labels=c(" < 100 "," < 1,000 "," < 10,0000 "," < 100,000 "," \u2265 100,000 "), 
        showNA=F
    ) +
    tm_shape(cshapes_proj_non_eu) + 
    tm_polygons(col = "white", lwd=.35) +
    tm_shape(cshapes_proj_euefta) + 
    tm_polygons(col = "grey",lwd=.35) +
    tm_layout(frame = FALSE,
        legend.position = c("right","bottom"),
        legend.text.size = 0.65,
        legend.width = 10,
        legend.title.fontfamily = "lm_regular",
        legend.text.fontfamily ="lm_math"
    )


# print map
print(map_N)

